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Abstract. - We study the dynamics of two symmetrically coupled populations of identical leaky 
integrate-and-fire neurons characterized by an excitatory coupling. Upon varying the coupling 
strength, we find symmetry-breaking transitions that lead to the onset of various chimera states 
as well as to a new regime, where the two populations are characterized by a different degree 
of synchronization. Symmetric collective states of increasing dynamical complexity are also ob- 
served. The computation of the the finite-amplitude Lyapunov exponent allows us to establish the 
chaoticity of the (collective) dynamics in a finite region of the phase plane. The further numeri- 
cal study of the standard Lyapunov spectrum reveals the presence of several positive exponents, 
indicating that the microscopic dynamics is high-dimensional. 
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Introduction. — Understanding the collective motion 
of ensembles / networks of oscillators is crucial in many con- 
texts, starting from neuronal circuits [T]. So far, most of 
the efforts have been devoted to the characterization of 
strong forms of synchronization. However, more subtle 
phenomena, like the onset of collective motion in an en- 
semble of (chaotic) units, which behave in a seemingly un- 
corrected way can also play a relevant role for information 
encoding. Collective chaos, meant as irregular dynamics 
of coarse-grained obscrvables, has been found in ensem- 
bles of fully coupled one-dimensional maps [2j[3] as well as 
in two-dimensional continuous-time oscillators [4-6]. In 
both classes of models, the single dynamical unit can be- 
have chaotically under the action of a periodic forcing (in 
non-invertible maps, there is even no need of a periodic 
forcing). What does it happen in ensembles of phase- 
oscillators which cannot become chaotic under the action 
of any forcing? The evolution of a (formally infinite) popu- 
lation of oscillators is ruled by a self-consistent (nonlinear) 
functional equation for the probability density. Given the 
infinite dimensionality of the model, the population could, 
in principle, behave chaotically, irrespective of the "struc- 
ture" of the single oscillators. In spite of this potentiality, 
only a few examples of low-dimensional chaotic collective 
motion have been found in ensembles of phase-oscillators 



[7]|8] . One reason is that most of the models so far investi- 
gated are based on sinusoidal force fields (in the following 
we refer to them as to sinusoidal oscillators); in this setup, 
there is little space for a high dimensional dynamics, since 
no matter how many oscillators are involved, there are 
always N — 3 constants of motion [8]. This "degener- 
acy" is not present in typical pulse-coupled networks of 
neurons, where different force fields arc usually assumed. 
A prototypical example is that of leaky integrate-and-fire 
(LIF) neurons, characterized by a linear force field. It 
is, in fact, not suprising that the first instance of a non- 
trivial collective motion has been found in an ensemble of 
LIF neurons. We refer to Partial Synchronization (PS) 
[5], a regime characterized by a periodic macroscopic dy- 
namics and a quasiperiodic microscopic motion, with the 
additional subtlety that the average inter-spike interval of 
the single neurons differs from the period of the collective 
variable. More recently, this type of behaviour has been 
observed also in a population of sinusoidal oscillators, in 
the presence of a suitable nonlinear coupling |10j . 

The only and quite striking evidence of an irregular col- 
lective dynamics has been recently found in an ensemble 
of LIF neurons in the presence of a random distribution of 
the input currents [TT] (this setup, where the single neu- 
rons are characterized by different spiking rates, is anal- 
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ogous to that of the usual Kuramoto model, where the 
single oscillators have different bare frequencies). On the 
other hand, the model studied in [TT] has the further pe- 
culiarity of exhibiting a negative maximal Lyapunov ex- 
ponent it is, in fact, an example of stable chaos |12j . 
However, in the absence of disorder, no example of irreg- 
ular dynamics has yet been found. 

A slightly more complex but meaningful setup is that 
of two symmetrically coupled populations of identical os- 
cillators. This is the simplest instance of "network-of- 
networks" that is often invoked as a paradigm for neural 
systems [T3]. With reference to sinusoidal oscillators, this 
setup has revealed the onset of chimera states (one of the 
two populations is fully synchronized, while the oscillators 
of the other one are not synchronized at all [Hj), as well 
as more complex macroscopic states with periodic |15| and 
quasi-periodic [TB] collective oscillations. In the present 
Letter we study the two-population setup with reference to 
LIF neurons for different values of the coupling strengths 
between and within the two populations. We find various 
kinds of symmetry broken states some of which are simi- 
lar to those observed in JTSHU] and a new one, where the 
two populations are both partially synchronized, but with 
a different degree. More interesting is the parameter re- 
gion where the collective motion is chaotic, as indicated by 
the finite-amplitude Lyapunov exponent (FALE) [17] and 
confirmed by the computation of the standard Lyapunov 
spectrum which reveals the existence of several positive 
exponents. 

The model. — We consider two fully coupled net- 
works, each made of N LIF oscillators. Following 
Refs. [T5], the membrane potential Xj(t) of the j — th 
oscillator (j = 1, . . . , N) of the fcth population (k = 0, 1) 
evolves according to the differential equation, 



,(fc) 
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xf } (t)+g s E^(t)+g c E^- k Ht) (1) 



where a > 1 is the suprathrcshold input current, while 
g s > and g c > gauge the self- and, resp., cross- 
coupling strength of the excitatory interaction. Whenever 
the membrane potential reaches the threshold — 1, 

(k) 

it is reset to =0, while a so-called cv-pulse is sent 
and instantaneously received by all the neurons. The field 
E^ k \t) represents the linear superposition of the pulses 
emitted within the population k in the past. It can be 
shown [T5] that E^ (£) satisfies the differential equation 



E (k \t) + 2aE^(t) + a 2 E^(t) 
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where tj „ is the nth spiking time of the jth neuron within 
the population k, and the sum is restricted to times smaller 
than t. In the limit case g s = g c = g, the two populations 
can be seen as a single one made of 2N neurons with an 
effective coupling constant G = 2g. 

The degree of synchronization can be quantified by 
introducing the typical order parameter used for phase 
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Fig. 1: (Color Online) Phase diagram in the (g c , gr s )-plane of 
the model (II 1 12 p for a = 1.3 and a = 9. FS indicates Full 
Synchronization (both populations fire at the same time); PS- 
FS indicates that the two populations are is the FS and PS 
regimes, respectively; PS1-PS2 indicates that both populations 
are in a PS regime, although with a different degree of synchro- 
nization; APS indicates Antiphase Partial Synchronization, i.e. 
the two fields exhibit the same behaviour though being in an- 
tiphase; TORUS indicates a collective quasi-periodic motion; 
finally, CHAOS indicates collective chaotic motion. 
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where 9^ is the 



phase of the jth oscillator, that can be properly defined 
by suitably rescaling the time variable [19], 9 A (t) — 



2ir(t— t) K n)/(tq^n — K where n identifies the last spike 



emitted by the jth neuron, while q indicates the neuron 
that has emitted the last spike. One can verify that this 
phase is bounded between and 2tt, as it should. It is in- 
teresting to see that the application of this definition to the 
PS regime described by van Vreeswijk [9] reveals that the 
order parameter fluctuates periodically. In other words PS 
differs from the regime observed in the Kuramoto model 
above the synchronization threshold, where the order pa- 
rameter is constant in time [20] , 

Phase Diagram. — The equations have been inte- 
grated by extending the event-driven approach described 
e.g. in [TB]. In practice, the (linear) equations of motion 
are solved analytically in between two consecutive spike- 
emissions, obtaining a suitable map. Since the ordering of 
the single-neuron potentials does not change within each 
population, the next firing event can be easily determined 
by comparing the neurons that are closest to threshold 
within each of two populations. In spite of the concep- 
tual simplicity and the effectiveness of the code, one must 
be nevertheless careful in handling nearly singular cases, 
when many neurons almost cluster together. In order to 
avoid the spurious clustering, due to numerical roundoff, 
we have changed variables, introducing and monitoring 
the logarithm of the difference of the membrane poten- 
tials of two successive neurons. This requires some care 
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Fig. 2: The macroscopic attractors displayed by reporting the 
fields £ (0) vs£ (1) for four different non chaotic phases, namely 
(a) PS-FS (<?c = 0.07, g a = 0.1), (b) PS1-PS2 (g c = 0.02, g a = 
0.17), (c) APS (g c = 0.07, g a = 0.35) and (d) TORUS (g c = 

0. 07, g 3 = 0.3) (for the exact definitions see the text). The 
grey curves reported in the panel (d) are the Poincare sections 
obtained by imposing that the sum 

£(0) _|_ i s maximal. 

in defining the right number of variables: given N poten- 
tials, one naturally has TV — 1 differences that have to be 
complemented by a proper TVth variable (for more details 
see [21]). 

The phase plane (g c ,9s) shown in Fig. [T]has been ob- 
tained by studying the model (|1I2|) for a = 1.3 and a = 9. 
The diagram is semiquantitative in the sense that a much 
more detailed work would be needed to identify exactly 
the stability borders of the different regimes. Along the 
diagonal (<? = g s = g c ) the model reduces to that for a 
single population with coupling strength G = 2g. For our 
choice of a and a values, the system exhibits PS, since we 
are below the critical value Go = 0.425 [18] above which 
the splay state is stable (the splay state is a regime charac- 
terized by a constant spiking rate and thereby a constant 
field, i.e. no collective dynamics). Below the diagonal, the 
evolution is still symmetric but fully synchronized (FS), 

1. e. all neurons of both populations fire together. More 
intriguing is the region above the diagonal, that is charac- 
terized by a spontaneous symmetry breaking: one popula- 
tion fully synchronizes, while the other is in a PS regime, 
i.e. we are in presence of a generalized chimera state (here 
termed PS-FS). This can be appreciated by looking at the 
synchronization parameter r^ k ' of the two populations, one 
of which is equal to one, while the other oscillates periodi- 
cally close to 0.8. By following [15], this state can be classi- 
fied as a periodically breathing chimera. In this regime, the 
two populations are characterized by a microscopically pe- 
riodic and quasi-periodic behaviour, respectively. In spite 
of this qualitative difference, the two (macroscopic) fields 
EM and E<>V are both periodic and phase locked (see 
Fig. [2k). This means that the neurons subject to two 
different linear combinations of E^ and E^> behave dif- 



ferently: a population locks with the forcing field, while 
the other one behaves quasi-periodically. Another even 
more interesting symmetry broken state can be observed 
for larger ^-values and g c < 0.055; in this case both pop- 
ulations exhibit PS, but their dynamics take place over 
two different attractors with two different degrees of syn- 
chronization (PS1-PS2 regime), as shown in Fig. [3k. Like 
in the PS-FS regime, the two fields behave periodically 
(with the same period) and are phase locked, as it can 
be appreciated by looking at the closed curve 

E (o) 

versus 

£?(!) in Fig. [5]3. However, at variance with PS-FS, here 
both populations exhibit quasi-periodic motions. In other 
words we are in presence of a different symmetry break- 
ing, where two populations with distinct quasi-periodic 
motions spontaneously emerge. For yet larger g s values, 
the equivalence between the collective dynamics of the two 
population is restored, the only difference being a phase 
shift between the two fields, which oscillate in antiphase 
and this is why we term this regime Antiphase Partial 
Synchronization (APS). In the APS phase, for finite N, 
the instantaneous maximum Lyapunov exponent strongly 
fluctuates and we cannot rule out the possible existence of 
some form of weak chaos, analogous to the one discussed 
in [52] for a model of diluted neural network, i.e. a chaotic 
behaviour that disappears in the thermodynamic limit. In 
a strip above the chaotic region (discussed below), one can 
observe collective quasipcriodic motion. This means that 
the quasipcriodic motion of the fields is accompanied by 
a dynamics of the single neurons along a torus T 3 . An 
analogous regime has been previously reported in [5] in 
the context of a population of coupled Stuart-Landau os- 
cillators. Here, we find it in a model where the single units 
are described by a single variable. Furthermore, we have 
characterized the motion on the macroscopic T 2 attrac- 
tor reported in Fig. [5] by estimating the winding num- 
bers for various values of the coupling strengths. We find 
that the winding number is typically quite small - on the 
order of ~ 10~ 2 but independent of the system size, 
indicating that the torus survives in the thermodynamic 
limit. Finally, for yet larger p s -values both populations 
converge towards a splay state. This is not surprising, as 
wc already know that for the chosen a- and a-values, the 
splay state is stable in a single population of neurons for 
G > G = 0.425. 

Collective Chaos. — In a limited region above the 
diagonal and for g c > 0.055 the collective behaviour is 
irregular, this is clearly seen by observing the two macro- 
scopic attractors, corresponding to the two populations, 
reported in Fig. [3]d. Nonetheless, from Fig. [3J it ap- 
pears that the associated order parameters behave almost 
periodically, but this periodicity is only apparent since it 
reflects only the larger time scale present in the system, 
the one associated to the modulation of the fields E^ k \ 
while the irregularity of the dynamics are more evident 
at smaller time scales. In particular, at least two other 
time scales are present: a scale 0(1) associated to the 
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Fig. 3: (Color Online) Macroscopic attractors displayed by re- 
porting P = E + aE vs E for a PS1-PS2 state (a) and a 
chaotic phase (b), the time evolution of the corresponding or- 
der parameters r 1 - ' and r' 1 ^ is also reported in (c) and (d). The 
variables corresponding to population (resp. 1) are shown as 
black solid lines (resp. red dashed lines) in (c) and (d); while 
in (a) the internal black (resp. external red) curve refers to 
population (resp. 1) and in (b) an unique attractor has been 
reported for clarity reasons, since the two attractors are over- 
lapping. The data reported in (a) and (c) refer to g c — 0.02 
and g 3 — 0.17, while those shown in (b),(d) to g c — 0.08 and 
g s = 0.16. 
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Fig. 4: Finite amplitude Lyapunov exponents Af versus the 
logarithm of the perturbation amplitude A for three differ- 
ent system sizes: namely, N — 400 (black circles), N = 800 
(green squares) and N = 1,600 (blue triangles). The am- 
plitudes A have been estimated by considering the euclidean 
distance among the perturbed and unperturbed fields. The 
dashed (magenta) line indicates the maximal microscopic Lya- 
punov exponent Ai for N = 1,600 obtained by following the 
orbit over a time span containing 10 spikes and after discard- 
ing a transient composed by 10 6 spikes. Af has been estimated 
by averaging over 25, 000 — 50, 000 different trajectories. The 
results have been obtained for g c — 0.08 and g s = 0.16. 



firing period of each specific neuron and a scale 0(1/N) 
corresponding to the interspike interval between two suc- 
cessive spike emissions in the network. In order to make 
a quantitative assessment of the chaoticity we have first 
studied the FALE Af- The FALE can be determined from 
the growth rate of a small finite perturbation for different 
amplitudes A of the perturbation itself (after averaging 
over different trajectories) [17] . This is done by randomly 
perturbing the coordinates (both fields and the membrane 
potentials of the two populations) of a generic configuara- 
tion on the attractor. The results for g s = 0.16, g c = 0.08 
and three different system sizes are plotted in Fig. For 
small A values, Af grows with A, since the perturbation 
needs first to converge towards the most expanding direc- 
tion, while the final drop is a manifestation of the satu- 
ration of the perturbation amplitude. It is the height of 
the intermediate plateau which measures the amplitude of 
the FALE. Since the height is independent of N, one can 
conjecture that the collective motion is chaotic and stays 
chaotic in the thermodynamic limit. Besides the data re- 
ported in Fig. SI we have checked that the plateau height 
is not influenced by changing the amplitude of the initial 
perturbation and by performing tests up to iV=6,400. 

It is instructive to compare A^ with the maximum Ai of 
the standard Lyapunov spectrum. For small A values, the 
two indicators should coincide, but there is no reason for 
the agreement to persist at larger amplitudes. In fact, a 
second plateau has been detected in the context of globally 
coupled maps [5J[3]. The plateau occurring for larger A's 
has been interpreted as an indication that the chaoticity 
of collective variables differs from that of the microscopic 
ones. Since, in the present case, we observe only a single 
plateau, it is crucial to verify its compatibility with Ai. 
The horizontal dashed line in Fig. 2] corresponds to the 
Lyapunov exponent determined for A=l,600 (for the de- 
pendence of Ai on N, see below). The two indicators are 
consistent with each other and this means that collective 
variables are as chaotic as the microscopic ones. 

Having established the existence of one unstable direc- 
tion, the next question is to determine how many such di- 
rections are present. Unfortunately, the concept of FALE 
allows to determine just one exponent. As a consequence, 
we turn our attention to the usual Lyapunov spectrum, 
well aware that the "microscopic" exponents do not neces- 
sarily reproduce the chaoticity of the "macroscopic" vari- 
ables. 

In the present model, the Lyapunov spectrum {Ai} is 
composed of N + 3 exponents (i = 1, . . . , N + 3). The 
phase space dimension is in fact equal to N + 4 (N po- 
tentials plus 2 two-dimensional equations for the fields), 
but the direction corresponding to the zero-Lyapunov ex- 
ponent is implicitly eliminated as a consequence of taking 
the Poincare section [18] . In Fig. [5] we have plotted the 
first part of the spectrum (with the exception of Ai, for the 
clarity of presentation) with the usual normalization i/N 
of the x axis. The figure clearly indicates that the spec- 
trum becomes increasingly flat and converges to zero. This 
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Fig. 5: (Color Online) Lyapunov spectra Ai versus i/N for 
three N- values: namely, N — 50 (filled black circles), N = fOO 
(empty red squares) and N = 200 (empty green triangles). In 
the inset the first three Lyapunov exponents are reported as a 
function of N: Ai (blue circles), A2 (turquoise squares) and A3 
(magenta triangles). The Lyapunov esponents have been ob- 
tained by following the dynamics in the real and tangent space 
for a time span containing f 8 — 5 x 10 9 spikes, after discarding 



a transient period of f 6 spikes. 
g c = 0.08 and g 3 = 0.16. 
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can be understood in the following way. In the thermo- 
dynamic limit, the dynamics of globally coupled identical 
oscillators can be viewed as that of single oscillators forced 
by the same (self-consistent) field [S]. As a result, in a first 
approximation, we expect all Lyapunov exponents to be 
equal to the conditional Lypunov exponent A c obtained 
by forcing a single LIF neuron (identical to the others) 
with the self-consistent fields (obtained by integrating the 
whole ensemble). We have found that A c = (within nu- 
merical accuracy). This justifies why the increasingly flat 
Lyapunov spectrum converges towards zero. Moreover, 
since a LIF neuron is described by a single variable, under 
no circumstance, A c can be strictly positive. On the other 
hand, A c can be negative, and we expect this to occur 
whenever the given population exhibits full synchroniza- 
tion. Since we know that no synchronization is observed 
in the chaotic regime, we can conclude that A c is not just 
very small, but it must be exactly equal to zero. 

Having established that the Lyapunov spectrum con- 
verges to zero, it is interesting to investigate its scaling be- 
haviour. By comparing the spectra obtained for different 
system sizes, it turns out that they scale as 1 /A* 3 , although 
we are unable to extrapolate the value of /? from the analy- 
sis of the relatively small systems that we have simulated. 
By comparing the spectra obtained for N = 100 and 200 
we can at most guess that f3 rj 1.5. This value is not far 
from (3 = 2, found analytically while studying the splay- 
state stability in single populations of LIF neurons [53] 
and numerically for the PS state [2"2"] . 

Finally, we turn our attention to the first part of the 
spectrum, where the flatness hypothesis does not hold. 



More precisely, we investigate the iV-dependence of the 
first three Lyapunov exponents which can be computed 
for larger lattices (up to N= 1,600). The maximal Lya- 
punov exponent appears to converge to a finite asymptotic 
value Ai = 0.0195(3). On the other hand, the second and 
third exponents grow systematically with N, both becom- 
ing positive for N > 200, with no clear evidence of an 
eventual saturation. These results suggest that the mi- 
croscopic chaos is high-dimensional (there is no reason to 
believe that the number of positive exponents is just equal 
to three). However, we cannot tell whether the number of 
positive exponents is extensive (proportional to N) or sub- 
extensive. 

Conclusions. — We have studied two symmetrically 
coupled populations of leaky integrate- and- fire neurons for 
different values of the self- (g s ) and cross-(g c ) coupling- 
strength. Some of the collective phenomena that we have 
identified are quite similar to those observed in the two- 
population setup of Kuramoto-like oscillators [53] . This is 
not surprising, since it is known that an ensemble of LIF 
neurons is equivalent, in the weak coupling limit, to the 
Kuramoto model [22], the only difference being that the 
coupling function is not purely sinusoidal. The onset of 
PS in both classes of models suggests that the equivalence 
can be extended to larger coupling strengths. However, 
since PS can be obtained in the Kuramoto setup only by 
invoking a more general kind of coupling |10j . it is legiti- 
mate to conclude that the relationship is more complicated 
than that suggested by the study of the weak-coupling 
limit. In fact, in this Letter we have found new dynami- 
cal regimes, such as a different PS dynamics for the two 
populations. A yet more intriguing phenomenon is the 
collective chaos that we have recognized as such from the 
computation of the finite-amplitude Lyapunov exponent. 
Altogether, a general question still stands: To what ex- 
tent are pulse-coupled oscillators equivalent to Kuramoto- 
like models? The identification of the mutual relationship 
would be highly beneficial in both areas. 

Another still open question is that of the degree of 
chaoticity of the collective dynamics. This problem is also 
connected to that of the asymptotic structure of the Lya- 
punov spectrum in the thermodynamic limit. A rough 
argument suggests that the spectrum should be flat and 
this is indeed approximately seen in the numerics. How- 
ever, the evident deviations observed in the vicinity of 
the maximum strongly suggest that the argument need be 
refined. Altogether, we observe that the Lyapunov spec- 
trum includes the FALE. This means that the evolution of 
not-so-small perturbations does not add anything, to that 
of infinitesimal ones and implies that the standard Lya- 
punov analysis is rich enough to account for the collective 
behaviour as well. This was not a priori obvious. The 
computation of the FALE Xp in globally coupled maps [3J 
reveals a different scenario, where Xf differs from the max- 
imum Lyapunov exponent. On the other hand, the more 
recent study of globally coupled Stuart-Landau oscillators 
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[26] provides an example where, like here, the chaoticity of 
the collective motion can be inferred from the Lyapunov 
spectrum. Moreover, what can we say about the role of 
the second and third exponents that are found to be pos- 
itive as well? Do they contribute to the microscopic dy- 
namics only, or also to the macroscopic one? A careful 
analysis based on the study of the corresponding covari- 
ant Lyapunov vectors |27] might help to clarify this point. 
An alternative and more direct approach could be that 
of (numerically) integrating the self-consistent dynamical 
equation for the probability densities of the membrane po- 
tentials in the two families. However, it is not easy to 
pursue this latter perspective: because of the occasional 
formation of strongly clusterized states, it is necessary to 
partition the phase space into a huge number of small cells. 
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